function [x,it,nrm,xx] = redblack2(ne,dc,sa,ns,sw,x,tol,maxkit)

global bcL bcR
n       = length(x);
j1      = x(1:end-1, 1); 
k1      = x(end, 1);
jr = zeros(ne*2,1);
jb = jr;
for i = 1:4:ne*2
   jr(i)      = 1;
   jr(i+1)    = 1;
end
for i = 1:4:((ne-mod(ne,2))*2)
   jb(i+2)    = 1;
   jb(i+3)   = 1;
end
I_r = diag(jr);
I_b = diag(jb);
k0      = 0;
it      = 1;
jwarn   = 0;
jjit    = 0;
maxjit  = 10000; 
nrm(1) = norm( feval('respfct',[j1' k1]') );
xx(:,1)=x;
while nrm(it) > tol(1)*nrm(it)+tol(2) && it <= maxkit  
    % get response matrices for current k value
    j0  = zeros(ne*2,1);
    [R F A L M] = respmtx(ne, dc, sa, ns, sw, [j1' k1]');
    R = sparse(R); M = sparse(M); MR = M*R;
    jit = 0;
    while norm(j1-j0) > tol(2) && jit < maxjit
        sumj    = sum(R*j1);
        j0      = j1;
        j1      = MR*j1/sumj;
        jit     = jit + 1;
    end
    if (jit==maxjit)
        jjit    = jjit+1;
        jwarn   = max(1,jwarn);
    end
    k0 = k1;
    % keep in mind j's are INCOMING currents
    jout    = R*j1;
    gain    = F*j1;
    loss    = A*j1 + jout(1)*(bcL==0) + jout(end)*(bcR==0);
    k1      = gain/loss;
    it      = it + 1; 
    nrm(it) = norm( feval('respfct',[jout' k1]') ); 
    xx(:,it)=[jout' k1]';
end

jout(:,1)=jout(:,1)/sum(jout(:,1));
if(jwarn>0)
disp('warning: at least one inner did not converge')
disp(['# of kits: ',num2str(it),' # of jjits ',num2str(jjit)])
end

x = [jout' k1]'; % the final output

